1 Goal

Preparation of figures for publication

2 Approach

library(Seurat)
library(viridis)
library(org.Hs.eg.db)
library(AnnotationDbi)
library(tidyverse)
library(purrr)
library(patchwork)
library(grid)
library(clusterProfiler)
set.seed(42)
boost <- readRDS("../data/boost_data_seurat.rds")
cds <- readRDS("../data/boost_data_cds.rds")

3 Figures

cols.cluster <- hcl.colors(6, "Spectral") %>% rev()
cols.cluster <- cols.cluster[c(1, 3, 4, 5, 6, 2)]

3.1 Figure 1 F: Dimensional reduction plot in UMAP space, cells colored by cluster membership, or pseudotime

DimPlot(boost, cols = cols.cluster, pt.size = 2) + theme(axis.text = element_text(size = 25), 
    axis.title = element_text(size = 30), plot.title = element_blank(), legend.text = element_text(size = 20)) + 
    guides(colour = guide_legend(override.aes = list(size = 8))) + scale_y_continuous(breaks = c(-5, 
    0, 5))

FeaturePlot(boost, features = "pseudotime", cols = plasma(length(boost$pseudotime)), 
    pt.size = 2) + theme(axis.text = element_text(size = 25), axis.title = element_text(size = 30), 
    plot.title = element_blank(), legend.text = element_text(size = 20), legend.key.size = unit(0.8, 
        "cm")) + scale_y_continuous(breaks = c(-5, 0, 5))
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.

3.2 Figure 1 G: Heatmap of specific marker genes. Cells are ordered by pseudotime, with an additional colorbar for the respective cluster.

genes <- c("MKI67", "SOX2", "HES1", "HES5", "NES", "PAX6", "TUBB3", "STMN1", "FAT3", 
    "DCX", "RBFOX3", "MAPT", "SNAP25", "CAMK2A", "CAMK4", "NRXN1", "NCAM1", "DLG4", 
    "SYP", "SYN1", "SYT1", "ANK3", "SPTBN4", "SCN2A", "GRIN1", "GRIN2B", "GRIA1", 
    "GRIA2", "GABRA3", "GABRB2", "ARC", "NPAS4", "FOXG1")
DoMultiBarHeatmap2 <- function(object, features = NULL, cells = NULL, group.by = "ident", 
    additional.group.by = NULL, additional.group.sort.by = NULL, cols.use = NULL, 
    group.bar = TRUE, disp.min = -2.5, disp.max = NULL, slot = "scale.data", assay = NULL, 
    label = TRUE, size = 5.5, hjust = 0, angle = 45, raster = TRUE, draw.lines = TRUE, 
    lines.width = NULL, group.bar.height = 0.02, combine = TRUE, group.name = NULL, 
    additional.group.name = NULL) {
    cells <- cells %||% colnames(x = object)
    if (is.numeric(x = cells)) {
        cells <- colnames(x = object)[cells]
    }
    assay <- assay %||% DefaultAssay(object = object)
    DefaultAssay(object = object) <- assay
    features <- features %||% VariableFeatures(object = object)
    features <- rev(x = unique(x = features))
    disp.max <- disp.max %||% ifelse(test = slot == "scale.data", yes = 2.5, no = 6)
    possible.features <- rownames(x = GetAssayData(object = object, slot = slot))
    if (any(!features %in% possible.features)) {
        bad.features <- features[!features %in% possible.features]
        features <- features[features %in% possible.features]
        if (length(x = features) == 0) {
            stop("No requested features found in the ", slot, " slot for the ", assay, 
                " assay.")
        }
        warning("The following features were omitted as they were not found in the ", 
            slot, " slot for the ", assay, " assay: ", paste(bad.features, collapse = ", "))
    }
    
    if (!is.null(additional.group.sort.by)) {
        if (any(!additional.group.sort.by %in% additional.group.by)) {
            bad.sorts <- additional.group.sort.by[!additional.group.sort.by %in% 
                additional.group.by]
            additional.group.sort.by <- additional.group.sort.by[additional.group.sort.by %in% 
                additional.group.by]
            if (length(x = bad.sorts) > 0) {
                warning("The following additional sorts were omitted as they were not a subset of additional.group.by : ", 
                  paste(bad.sorts, collapse = ", "))
            }
        }
    }
    
    data <- as.data.frame(x = as.matrix(x = t(x = GetAssayData(object = object, slot = slot)[features, 
        cells, drop = FALSE])))
    
    object <- suppressMessages(expr = StashIdent(object = object, save.name = "ident"))
    group.by <- group.by %||% "ident"
    groups.use <- object[[c(group.by, additional.group.by[!additional.group.by %in% 
        group.by])]][cells, , drop = FALSE]
    plots <- list()
    for (i in group.by) {
        data.group <- data
        if (!is_null(additional.group.by)) {
            additional.group.use <- additional.group.by[additional.group.by != i]
            if (!is_null(additional.group.sort.by)) {
                additional.sort.use = additional.group.sort.by[additional.group.sort.by != 
                  i]
            } else {
                additional.sort.use = NULL
            }
        } else {
            additional.group.use = NULL
            additional.sort.use = NULL
        }
        
        group.use <- groups.use[, c(i, additional.group.use), drop = FALSE]
        
        for (colname in colnames(group.use)) {
            if (!is.factor(x = group.use[[colname]])) {
                group.use[[colname]] <- factor(x = group.use[[colname]])
            }
        }
        
        if (draw.lines) {
            lines.width <- lines.width %||% ceiling(x = nrow(x = data.group) * 0.0025)
            placeholder.cells <- sapply(X = 1:(length(x = levels(x = group.use[[i]])) * 
                lines.width), FUN = function(x) {
                return(Seurat:::RandomName(length = 20))
            })
            placeholder.groups <- data.frame(rep(x = levels(x = group.use[[i]]), 
                times = lines.width))
            group.levels <- list()
            group.levels[[i]] = levels(x = group.use[[i]])
            for (j in additional.group.use) {
                group.levels[[j]] <- levels(x = group.use[[j]])
                placeholder.groups[[j]] = NA
            }
            
            colnames(placeholder.groups) <- colnames(group.use)
            rownames(placeholder.groups) <- placeholder.cells
            
            group.use <- sapply(group.use, as.vector)
            rownames(x = group.use) <- cells
            
            group.use <- rbind(group.use, placeholder.groups)
            
            for (j in names(group.levels)) {
                group.use[[j]] <- factor(x = group.use[[j]], levels = group.levels[[j]])
            }
            
            na.data.group <- matrix(data = NA, nrow = length(x = placeholder.cells), 
                ncol = ncol(x = data.group), dimnames = list(placeholder.cells, colnames(x = data.group)))
            data.group <- rbind(data.group, na.data.group)
        }
        
        order_expr <- paste0("order(", paste(c(i, additional.sort.use), collapse = ","), 
            ")")
        group.use = with(group.use, group.use[eval(parse(text = order_expr)), , drop = F])
        
        plot <- Seurat:::SingleRasterMap(data = data.group, raster = raster, disp.min = disp.min, 
            disp.max = disp.max, feature.order = features, cell.order = rownames(x = group.use), 
            group.by = group.use[[i]])
        
        if (group.bar) {
            pbuild <- ggplot_build(plot = plot)
            group.use2 <- group.use
            cols <- list()
            na.group <- Seurat:::RandomName(length = 20)
            for (colname in rev(x = colnames(group.use2))) {
                if (colname == i) {
                  colid = group.name
                } else {
                  colid = additional.group.name
                }
                
                # Default
                cols[[colname]] <- c((scales::hue_pal())(length(x = levels(x = group.use[[colname]]))))
                
                # Overwrite if better value is provided
                if (!is_null(cols.use[[colname]])) {
                  req_length = length(x = levels(group.use))
                  if (length(cols.use[[colname]]) < req_length) {
                    warning("Cannot use provided colors for ", colname, " since there aren't enough colors.")
                  } else {
                    if (!is_null(names(cols.use[[colname]]))) {
                      if (all(levels(group.use[[colname]]) %in% names(cols.use[[colname]]))) {
                        cols[[colname]] <- as.vector(cols.use[[colname]][levels(group.use[[colname]])])
                      } else {
                        warning("Cannot use provided colors for ", colname, " since all levels (", 
                          paste(levels(group.use[[colname]]), collapse = ","), ") are not represented.")
                      }
                    } else {
                      cols[[colname]] <- as.vector(cols.use[[colname]])[c(1:length(x = levels(x = group.use[[colname]])))]
                    }
                  }
                }
                
                # Add white if there's lines
                if (draw.lines) {
                  levels(x = group.use2[[colname]]) <- c(levels(x = group.use2[[colname]]), 
                    na.group)
                  group.use2[placeholder.cells, colname] <- na.group
                  cols[[colname]] <- c(cols[[colname]], "#FFFFFF")
                }
                names(x = cols[[colname]]) <- levels(x = group.use2[[colname]])
                
                y.range <- diff(x = pbuild$layout$panel_params[[1]]$y.range)
                y.pos <- max(pbuild$layout$panel_params[[1]]$y.range) + y.range * 
                  0.015
                y.max <- y.pos + group.bar.height * y.range
                pbuild$layout$panel_params[[1]]$y.range <- c(pbuild$layout$panel_params[[1]]$y.range[1], 
                  y.max)
                
                plot <- suppressMessages(plot + annotation_raster(raster = t(x = cols[[colname]][group.use2[[colname]]]), 
                  xmin = -Inf, xmax = Inf, ymin = y.pos, ymax = y.max) + annotation_custom(grob = grid::textGrob(label = colid, 
                  hjust = 0, gp = gpar(cex = 0.75)), ymin = mean(c(y.pos, y.max)), 
                  ymax = mean(c(y.pos, y.max)), xmin = Inf, xmax = Inf) + coord_cartesian(ylim = c(0, 
                  y.max), clip = "off"))
                
                if ((colname == i) && label) {
                  x.max <- max(pbuild$layout$panel_params[[1]]$x.range)
                  x.divs <- pbuild$layout$panel_params[[1]]$x.major %||% pbuild$layout$panel_params[[1]]$x$break_positions()
                  group.use$x <- x.divs
                  label.x.pos <- tapply(X = group.use$x, INDEX = group.use[[colname]], 
                    FUN = median) * x.max
                  label.x.pos <- data.frame(group = names(x = label.x.pos), label.x.pos)
                  plot <- plot + geom_text(stat = "identity", data = label.x.pos, 
                    aes_string(label = "group", x = "label.x.pos"), y = y.max + y.max * 
                      0.03 * 0.5, angle = angle, hjust = hjust, size = size)
                  plot <- suppressMessages(plot + coord_cartesian(ylim = c(0, y.max + 
                    y.max * 0.002 * max(nchar(x = levels(x = group.use[[colname]]))) * 
                      size), clip = "off"))
                }
            }
        }
        plot <- plot + theme(line = element_blank())
        plots[[i]] <- plot
    }
    if (combine) {
        plots <- CombinePlots(plots = plots)
    }
    return(plots)
}
boost_subset <- subset(boost, new_order %in% c("1", "2", "3", "4", "5"))
cols.use <- list(pseudotime = plasma(length(boost_subset$pseudotime)), ident = cols.cluster)
DoMultiBarHeatmap2(boost_subset, features = genes, group.by = "pseudotime", additional.group.by = "ident", 
    draw.lines = FALSE, label = FALSE, cols.use = cols.use, raster = FALSE) + scale_fill_viridis_c(option = "viridis") + 
    NoLegend() + theme(axis.text = element_text(size = 14, colour = "black")) + scale_y_discrete(position = "left") + 
    geom_hline(yintercept = c(3.5, 9.5, 12.5, 16.5, 23.5, 27.5), color = "white", 
        size = 1)
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
# DoHeatmap(boost_subset, features = genes, group.by='ident', draw.lines = FALSE,
# label=FALSE, raster = FALSE)+ scale_fill_viridis_c(option='viridis')+
# theme(legend.text = element_text(size = 20), legend.key.size = unit(0.8, 'cm'),
# legend.text.align = 1, legend.title = element_text(size=30))

3.3 Figure 1 H: NMI scores per cluster of neuronal cells.

sub <- subset(boost, idents = c("2", "3", "4", "5"))
a <- VlnPlot(sub, features = "NMIoverall", pt.size = 0.1, cols = cols.cluster[2:5]) + 
    NoLegend() + theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20), 
    plot.title = element_text(size = 18)) + xlab("") + ggtitle("NMI") + ylab("")
b <- VlnPlot(sub, features = "NMIdiscriminable", pt.size = 0.1, cols = cols.cluster[2:5]) + 
    NoLegend() + theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20), 
    plot.title = element_text(size = 18)) + xlab("") + ggtitle("Discriminating NMI") + 
    ylab("")
c <- VlnPlot(sub, features = "NMIactivity", pt.size = 0.1, cols = cols.cluster[2:5]) + 
    NoLegend() + theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20), 
    plot.title = element_text(size = 18)) + xlab("") + ggtitle("Functionality index") + 
    ylab("")
# pw <- (a/b/c)&theme(axis.text.x = element_text(angle = 0, hjust =
# 0.5),axis.text=element_text(size=20),
# plot.title=element_text(size=25))&xlab('')
# pw&theme(axis.text.x = element_text(angle = 0, hjust =
# 0.5),axis.text=element_text(size=20),
# plot.title=element_text(size=22))&xlab('')&ylab('')
b + scale_y_continuous(breaks = c(0.49, 0.52, 0.55))
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.

c

3.4 Figure 2 D: Violin Plots of CASP9, CASP3, CASP7

genes <- c("CASP7", "CASP9", "CASP3")
p <- VlnPlot(boost_subset, features = genes, group.by = "new_order", pt.size = 0.1, 
    cols = cols.cluster, combine = F)  #+NoLegend()+xlab('')#+theme(axis.text.x = element_text(angle = 0,  hjust = 0.5))
# patchwork <- wrap_plots(p, byrow = T)

3.4.1 CASP9

p[[2]] + scale_y_continuous(breaks = c(0, 0.3, 0.6, 0.9)) + theme(axis.text.x = element_text(angle = 0, 
    hjust = 0.5), axis.text = element_text(size = 20), plot.title = element_blank()) + 
    xlab("") + ylab("") + NoLegend()
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.

3.4.2 CASP3

p[[3]] + scale_y_continuous(breaks = c(0, 1, 2, 3), labels = c("0.0", "1.0", "2.0", 
    "3.0")) + theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20), 
    plot.title = element_text(size = 18)) + xlab("") + ylab("") + NoLegend()
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.

3.4.3 CASP7

p[[1]] + # scale_y_continuous(breaks=c(0.0, 0.3,0.6, 0.9))+
theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20), 
    plot.title = element_blank()) + xlab("") + ylab("") + NoLegend()

3.5 Figure 2 I: Violin Plots BAX, BCL2, BAK1

genes <- c("BAX", "BCL2", "BAK1")
p <- VlnPlot(boost_subset, features = genes, group.by = "new_order", pt.size = 0.1, 
    cols = cols.cluster)  #+NoLegend()+xlab('')#+theme(axis.text.x = element_text(angle = 0,  hjust = 0.5))

3.5.1 BAX

p[[1]] + scale_y_continuous(breaks = c(0, 1, 2), labels = c("0.0", "1.0", "2.0")) + 
    theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20), 
        plot.title = element_text(size = 18)) + xlab("") + ylab("") + NoLegend()
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.

3.5.2 BCL2

p[[2]] + scale_y_continuous(breaks = c(0, 0.4, 0.8, 1.2)) + theme(axis.text.x = element_text(angle = 0, 
    hjust = 0.5), axis.text = element_text(size = 20), plot.title = element_text(size = 18)) + 
    xlab("") + ylab("") + NoLegend()
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.

3.5.3 BAK1

p[[3]] + # scale_y_continuous(breaks=c(0.0, 1.0, 2.0), labels=c('0.0', '1.0', '2.0'))+
theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20), 
    plot.title = element_text(size = 18)) + xlab("") + ylab("") + NoLegend()

3.6 Figure 2 M: Violin Plot CYCS

genes <- "CYCS"
p <- VlnPlot(boost_subset, features = genes, group.by = "new_order", pt.size = 0.1, 
    cols = cols.cluster)  #+NoLegend()+xlab('')#+theme(axis.text.x = element_text(angle = 0,  hjust = 0.5))
p[[1]] + scale_y_continuous(breaks = c(0, 1, 2, 3), labels = c("0.0", "1.0", "2.0", 
    "3.0")) + theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20), 
    plot.title = element_text(size = 18)) + xlab("") + ylab("") + NoLegend()
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.

3.7 Figure 3 C: Violin Plots BIRC2, BIRC5, AREL1, BIRC4 (XIAP)

genes <- c("BIRC2", "XIAP", "BIRC5", "AREL1")
p <- VlnPlot(boost_subset, features = genes, group.by = "new_order", pt.size = 0.1, 
    cols = cols.cluster, combine = F)  #+NoLegend()+xlab('')#+theme(axis.text.x = element_text(angle = 0,  hjust = 0.5))

3.7.1 BIRC2

p[[1]] + # scale_y_continuous(breaks=c(0.0, 1.0, 2.0, 3.0), labels=c('0.0', '1.0', '2.0',
# '3.0'))+
theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20), 
    plot.title = element_text(size = 18)) + xlab("") + ylab("") + NoLegend()

3.7.2 BIRC5

p[[3]] + scale_y_continuous(breaks = c(0, 1, 2), labels = c("0.0", "1.0", "2.0")) + 
    theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20), 
        plot.title = element_text(size = 18)) + xlab("") + ylab("") + NoLegend()
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.

3.7.3 AREL1

p[[4]] + # scale_y_continuous(breaks=c(0.0, 1.0, 2.0, 3.0), labels=c('0.0', '1.0', '2.0',
# '3.0'))+
theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20), 
    plot.title = element_text(size = 18)) + xlab("") + ylab("") + NoLegend()

3.7.4 XIAP/BIRC4

p[[2]] + # scale_y_continuous(breaks=c(0.0, 1.0, 2.0), labels=c('0.0', '1.0', '2.0'))+
theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20), 
    plot.title = element_text(size = 18)) + xlab("") + ylab("") + NoLegend()

3.8 Figure 3 F: Violin Plots HTRA2, SIVA1, SIAH1

genes <- c("HTRA2", "SIVA1", "SIAH1", "USP11")
p <- VlnPlot(boost_subset, features = genes, group.by = "new_order", pt.size = 0.1, 
    cols = cols.cluster)  #+NoLegend()+xlab('')#+theme(axis.text.x = element_text(angle = 0,  hjust = 0.5))

3.8.1 HTRA2

p[[1]] + # scale_y_continuous(breaks=c(0.0, 1.0, 2.0, 3.0), labels=c('0.0', '1.0', '2.0',
# '3.0'))+
theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20), 
    plot.title = element_text(size = 18)) + xlab("") + ylab("") + NoLegend()

3.8.2 SIVA1

p[[2]] + scale_y_continuous(breaks = c(0, 1, 2), labels = c("0.0", "1.0", "2.0")) + 
    theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20), 
        plot.title = element_text(size = 18)) + xlab("") + ylab("") + NoLegend()
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.

3.8.3 SIAH1

p[[3]] + scale_y_continuous(breaks = c(0, 1, 2), labels = c("0.0", "1.0", "2.0")) + 
    theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20), 
        plot.title = element_text(size = 18)) + xlab("") + ylab("") + NoLegend()
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.

3.9 Figure 3 G: Violin Plot USP11

p[[4]] + scale_y_continuous(breaks = c(0, 1, 2, 3), labels = c("0.0", "1.0", "2.0", 
    "3.0")) + theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20), 
    plot.title = element_text(size = 18)) + xlab("") + ylab("") + NoLegend()
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.

3.10 Supplementary Figure 1 H: Pseudotime trajectories calculated on the per cell expression values for specific genes.

my_plot_genes <- function(cds_subset, min_expr = NULL, cell_size = 0.75, nrow = NULL, 
    ncol = 1, panel_order = NULL, color_cells_by = "pseudotime", trend_formula = "~ splines::ns(pseudotime, df=3)", 
    label_by_short_name = TRUE, vertical_jitter = NULL, horizontal_jitter = NULL, 
    line.size = NULL) {
    assertthat::assert_that(methods::is(cds_subset, "cell_data_set"))
    tryCatch({
        pseudotime(cds_subset)
    }, error = function(x) {
        stop(paste("No pseudotime calculated. Must call order_cells first."))
    })
    colData(cds_subset)$pseudotime <- pseudotime(cds_subset)
    if (!is.null(min_expr)) {
        assertthat::assert_that(assertthat::is.number(min_expr))
    }
    assertthat::assert_that(assertthat::is.number(cell_size))
    if (!is.null(nrow)) {
        assertthat::assert_that(assertthat::is.count(nrow))
    }
    assertthat::assert_that(assertthat::is.count(ncol))
    assertthat::assert_that(is.logical(label_by_short_name))
    if (label_by_short_name) {
        assertthat::assert_that("gene_short_name" %in% names(rowData(cds_subset)), 
            msg = paste("When label_by_short_name = TRUE,", "rowData must have a column of gene", 
                "names called gene_short_name."))
    }
    assertthat::assert_that(color_cells_by %in% c("cluster", "partition") | color_cells_by %in% 
        names(colData(cds_subset)), msg = paste("color_cells_by must be a column in the", 
        "colData table."))
    if (!is.null(panel_order)) {
        if (label_by_short_name) {
            assertthat::assert_that(all(panel_order %in% rowData(cds_subset)$gene_short_name))
        } else {
            assertthat::assert_that(all(panel_order %in% row.names(rowData(cds_subset))))
        }
    }
    assertthat::assert_that(nrow(rowData(cds_subset)) <= 100, msg = paste("cds_subset has more than 100 genes -", 
        "pass only the subset of the CDS to be", "plotted."))
    assertthat::assert_that(methods::is(cds_subset, "cell_data_set"))
    assertthat::assert_that("pseudotime" %in% names(colData(cds_subset)), msg = paste("pseudotime must be a column in", 
        "colData. Please run order_cells", "before running", "plot_genes_in_pseudotime."))
    if (!is.null(min_expr)) {
        assertthat::assert_that(assertthat::is.number(min_expr))
    }
    assertthat::assert_that(assertthat::is.number(cell_size))
    assertthat::assert_that(!is.null(size_factors(cds_subset)))
    if (!is.null(nrow)) {
        assertthat::assert_that(assertthat::is.count(nrow))
    }
    assertthat::assert_that(assertthat::is.count(ncol))
    assertthat::assert_that(is.logical(label_by_short_name))
    if (label_by_short_name) {
        assertthat::assert_that("gene_short_name" %in% names(rowData(cds_subset)), 
            msg = paste("When label_by_short_name = TRUE,", "rowData must have a column of gene", 
                "names called gene_short_name."))
    }
    assertthat::assert_that(color_cells_by %in% c("cluster", "partition") | color_cells_by %in% 
        names(colData(cds_subset)), msg = paste("color_cells_by must be a column in the", 
        "colData table."))
    if (!is.null(panel_order)) {
        if (label_by_short_name) {
            assertthat::assert_that(all(panel_order %in% rowData(cds_subset)$gene_short_name))
        } else {
            assertthat::assert_that(all(panel_order %in% row.names(rowData(cds_subset))))
        }
    }
    assertthat::assert_that(nrow(rowData(cds_subset)) <= 100, msg = paste("cds_subset has more than 100 genes -", 
        "pass only the subset of the CDS to be", "plotted."))
    f_id <- NA
    Cell <- NA
    cds_subset = cds_subset[, is.finite(colData(cds_subset)$pseudotime)]
    
    
    cds_exprs <- SingleCellExperiment::counts(cds_subset)
    cds_exprs <- Matrix::t(Matrix::t(cds_exprs)/size_factors(cds_subset))
    cds_exprs <- reshape2::melt(round(as.matrix(cds_exprs)))
    if (is.null(min_expr)) {
        min_expr <- 0
    }
    colnames(cds_exprs) <- c("f_id", "Cell", "expression")
    cds_colData <- colData(cds_subset)
    cds_rowData <- rowData(cds_subset)
    
    cds_exprs <- as.data.frame(merge(cds_exprs, cds_rowData, by.x = "f_id", by.y = "row.names"))
    
    
    cds_exprs <- as.data.frame(merge(cds_exprs, cds_colData, by.x = "Cell", by.y = "row.names"))
    
    
    
    cds_exprs$adjusted_expression <- cds_exprs$expression
    if (label_by_short_name == TRUE) {
        if (is.null(cds_exprs$gene_short_name) == FALSE) {
            cds_exprs$feature_label <- as.character(cds_exprs$gene_short_name)
            cds_exprs$feature_label[is.na(cds_exprs$feature_label)] <- cds_exprs$f_id
        } else {
            cds_exprs$feature_label <- cds_exprs$f_id
        }
    } else {
        cds_exprs$feature_label <- cds_exprs$f_id
    }
    cds_exprs$f_id <- as.character(cds_exprs$f_id)
    cds_exprs$feature_label <- factor(cds_exprs$feature_label)
    new_data <- data.frame(pseudotime = colData(cds_subset)$pseudotime)
    model_tbl = fit_models(cds_subset, model_formula_str = trend_formula)
    model_expectation <- model_predictions(model_tbl, new_data = colData(cds_subset))
    colnames(model_expectation) <- colnames(cds_subset)
    expectation <- plyr::ddply(cds_exprs, plyr::.(f_id, Cell), function(x) {
        data.frame(expectation = model_expectation[x$f_id, x$Cell])
    })
    cds_exprs <- merge(cds_exprs, expectation)
    cds_exprs$expression[cds_exprs$expression < min_expr] <- min_expr
    cds_exprs$expectation[cds_exprs$expectation < min_expr] <- min_expr
    if (!is.null(panel_order)) {
        cds_exprs$feature_label <- factor(cds_exprs$feature_label, levels = panel_order)
    }
    q <- ggplot(aes(pseudotime, expression), data = cds_exprs)
    if (!is.null(color_cells_by)) {
        q <- q + geom_point(aes_string(color = color_cells_by), size = I(cell_size), 
            position = position_jitter(horizontal_jitter, vertical_jitter))
        if (class(colData(cds_subset)[, color_cells_by]) == "numeric") {
            q <- q + viridis::scale_color_viridis(option = "C")
        }
    } else {
        q <- q + geom_point(size = I(cell_size), position = position_jitter(horizontal_jitter, 
            vertical_jitter))
    }
    q <- q + geom_line(aes(x = pseudotime, y = expectation), data = cds_exprs, size = line.size)
    q <- q + scale_y_log10() + facet_wrap(~feature_label, nrow = nrow, ncol = ncol, 
        scales = "free_y")
    if (min_expr < 1) {
        q <- q + expand_limits(y = c(min_expr, 1))
    }
    q <- q + ylab("Expression")
    q <- q + xlab("pseudotime")
    # q <- q + monocle_theme_opts()
    q
}
genes <- c("SOX2", "NES", "TUBB3", "DCX", "SNAP25", "SYN1")
cds_subset <- cds[, SummarizedExperiment::colData(cds)$ident %in% c("1", "2", "3", 
    "4", "5")]
## Loading required package: monocle3
## Loading required package: SingleCellExperiment
## Loading required package: SummarizedExperiment
## Loading required package: MatrixGenerics
## Loading required package: matrixStats
## 
## Attaching package: 'matrixStats'
## The following object is masked from 'package:dplyr':
## 
##     count
## The following objects are masked from 'package:Biobase':
## 
##     anyMissing, rowMedians
## 
## Attaching package: 'MatrixGenerics'
## The following objects are masked from 'package:matrixStats':
## 
##     colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
##     colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
##     colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
##     colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
##     colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
##     colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
##     colWeightedMeans, colWeightedMedians, colWeightedSds,
##     colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
##     rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
##     rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
##     rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
##     rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
##     rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
##     rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
##     rowWeightedSds, rowWeightedVars
## The following object is masked from 'package:Biobase':
## 
##     rowMedians
## Loading required package: GenomicRanges
## Loading required package: GenomeInfoDb
## 
## Attaching package: 'SummarizedExperiment'
## The following object is masked from 'package:Seurat':
## 
##     Assays
## 
## Attaching package: 'monocle3'
## The following objects are masked from 'package:Biobase':
## 
##     exprs, fData, fData<-, pData, pData<-
lineage_cds <- cds_subset[rowData(cds_subset)$gene_short_name %in% "SOX2", ]
p <- list()
p[[1]] <- my_plot_genes(lineage_cds, color_cells_by = "new_order", cell_size = 4, 
    line.size = 2) + theme_grey(base_size = 40, ) + NoLegend() + scale_color_manual(values = c("#584B9F", 
    "#BAEEAE", "#FCDE85", "#ED820A", "#A71B4B")) + theme(axis.title = element_blank(), 
    axis.text.x = element_blank(), strip.background = element_rect(fill = "white", 
        colour = NA)) + scale_y_log10(breaks = c(0.01, 0.1, 1, 10), labels = c("0.01", 
    "0.1", "1", "10"))
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
lineage_cds <- cds_subset[rowData(cds_subset)$gene_short_name %in% "NES", ]

p[[2]] <- my_plot_genes(lineage_cds, color_cells_by = "new_order", cell_size = 4, 
    line.size = 2) + theme_grey(base_size = 40) + NoLegend() + scale_color_manual(values = c("#584B9F", 
    "#BAEEAE", "#FCDE85", "#ED820A", "#A71B4B")) + theme(axis.title.x = element_blank(), 
    axis.text.x = element_blank(), strip.background = element_rect(fill = "white", 
        colour = NA))
lineage_cds <- cds_subset[rowData(cds_subset)$gene_short_name %in% "TUBB3", ]

p[[3]] <- my_plot_genes(lineage_cds, color_cells_by = "new_order", cell_size = 4, 
    line.size = 2) + theme_grey(base_size = 40) + NoLegend() + scale_color_manual(values = c("#584B9F", 
    "#BAEEAE", "#FCDE85", "#ED820A", "#A71B4B")) + theme(axis.title = element_blank(), 
    strip.background = element_rect(fill = "white", colour = NA))
lineage_cds <- cds_subset[rowData(cds_subset)$gene_short_name %in% "DCX", ]

p[[4]] <- my_plot_genes(lineage_cds, color_cells_by = "new_order", cell_size = 4, 
    line.size = 2) + theme_grey(base_size = 40) + NoLegend() + scale_color_manual(values = c("#584B9F", 
    "#BAEEAE", "#FCDE85", "#ED820A", "#A71B4B")) + theme(axis.title = element_blank(), 
    axis.text.x = element_blank(), strip.background = element_rect(fill = "white", 
        colour = NA))
lineage_cds <- cds_subset[rowData(cds_subset)$gene_short_name %in% "SNAP25", ]

p[[5]] <- my_plot_genes(lineage_cds, color_cells_by = "new_order", cell_size = 4, 
    line.size = 2) + theme_grey(base_size = 40) + NoLegend() + scale_color_manual(values = c("#584B9F", 
    "#BAEEAE", "#FCDE85", "#ED820A", "#A71B4B")) + theme(axis.title = element_blank(), 
    axis.text.x = element_blank(), strip.background = element_rect(fill = "white", 
        colour = NA))
lineage_cds <- cds_subset[rowData(cds_subset)$gene_short_name %in% "SYN1", ]

p[[6]] <- my_plot_genes(lineage_cds, color_cells_by = "new_order", cell_size = 4, 
    line.size = 2) + theme_grey(base_size = 40) + NoLegend() + scale_color_manual(values = c("#584B9F", 
    "#BAEEAE", "#FCDE85", "#ED820A", "#A71B4B")) + theme(axis.title = element_blank(), 
    strip.background = element_rect(fill = "white", colour = NA))
wrap_plots(p[1:3], ncol = 1)

wrap_plots(p[4:6], ncol = 1)

3.11 Supplementary Figure 2 B: Violin Plots CASP8, CASP10, CASP2, CASP6

genes <- c("CASP8", "CASP10", "CASP2", "CASP6")
p <- VlnPlot(boost_subset, features = genes, group.by = "new_order", pt.size = 0.1, 
    cols = cols.cluster, combine = F)  #+NoLegend()+xlab('')#+theme(axis.text.x = element_text(angle = 0,  hjust = 0.5))
# patchwork <- wrap_plots(p, byrow = T)

3.11.1 CASP8

p[[1]] + # scale_y_continuous(breaks=c(0.0, 0.3,0.6, 0.9))+
theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20), 
    plot.title = element_blank()) + xlab("") + ylab("") + NoLegend()

3.11.2 CASP10

p[[2]] + # scale_y_continuous(breaks=c(0.0, 0.3,0.6, 0.9))+
theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20), 
    plot.title = element_blank()) + xlab("") + ylab("") + NoLegend()

3.11.3 CASP2

p[[3]] + # scale_y_continuous(breaks=c(0.0, 0.3, 0.6, 0.9))+
theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20), 
    plot.title = element_blank()) + xlab("") + ylab("") + NoLegend()

3.11.4 CASP6

p[[4]] + # scale_y_continuous(breaks=c(0.0, 1.0, 2.0, 3.0), labels=c('0.0', '1.0', '2.0',
# '3.0'))+
theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20), 
    plot.title = element_blank()) + xlab("") + ylab("") + NoLegend()

3.12 Supplementary Figure 2 C: Violin Plot ACIN1

genes <- c("ACIN1")
p <- VlnPlot(boost_subset, features = genes, group.by = "new_order", pt.size = 0.1, 
    cols = cols.cluster, combine = FALSE)
p[[1]] + scale_y_continuous(breaks = c(0, 1, 2), labels = c("0.0", "1.0", "2.0")) + 
    theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20), 
        plot.title = element_blank()) + xlab("") + ylab("") + NoLegend()
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.

3.13 Supplementary Figure 3 E: Violin Plot XAF1

genes <- c("XAF1")
p <- VlnPlot(boost_subset, features = genes, group.by = "new_order", pt.size = 0.1, 
    cols = cols.cluster)  #+NoLegend()+xlab('')#+theme(axis.text.x = element_text(angle = 0,  hjust = 0.5))
p[[1]] + scale_y_continuous(breaks = c(0, 0.4, 0.8)) + theme(axis.text.x = element_text(angle = 0, 
    hjust = 0.5), axis.text = element_text(size = 20), plot.title = element_blank()) + 
    xlab("") + ylab("") + NoLegend()
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.

3.14 Supplementary Figure 4 C: Glut/Gaba split violin Plots CASP3, CASP9, CASP6, BIRC2, XIAP, AREL1, BCL2, BAK1, BAX, USP11, SIVA1, CYCS

rm(list = setdiff(ls(), c("boost", "cols.cluster")))
SLC17A7 <- WhichCells(boost, expression = SLC17A7 > 0)
SLC17A6 <- WhichCells(boost, expression = SLC17A6 > 0)
TBR1 <- WhichCells(boost, expression = TBR1 > 0)
Glut <- union(SLC17A6, SLC17A7)
Glut <- union(Glut, TBR1)
GAD1 <- WhichCells(boost, expression = GAD1 > 0)
GAD2 <- WhichCells(boost, expression = GAD2 > 0)
GAB <- union(GAD1, GAD2)
both_exp <- FetchData(boost, vars = c("SLC17A7", "SLC17A6", "TBR1", "GAD1", "GAD2"), 
    cells = intersect(Glut, GAB))
both_exp <- data.frame(Glut = rowMeans(both_exp[, 1:3]), Gaba = rowMeans(both_exp[, 
    4:5]))
both_exp$cell_id <- rownames(both_exp)
both_exp <- both_exp %>% mutate(type = if_else(Glut > Gaba, "Glut", "Gaba"))
both_GAB <- filter(both_exp, type == "Gaba")
both_Glut <- filter(both_exp, type == "Glut")
GAB <- setdiff(GAB, Glut)
GAB <- c(GAB, both_GAB$cell_id)
Glut <- setdiff(Glut, GAB)
Glut <- c(Glut, both_Glut$cell_id)
boost[["old.ident"]] <- Idents(object = boost)
boost <- StashIdent(object = boost, save.name = "old.ident")
## With Seurat 3.X, stashing identity classes can be accomplished with the following:
## boost[["old.ident"]] <- Idents(object = boost)
Idents(object = boost, cells = GAB) <- "Gaba"
Idents(object = boost, cells = Glut) <- "Glut"
Rest <- WhichCells(boost, invert = TRUE, cells = c(GAB, Glut))
Idents(object = boost, cells = Rest) <- "rest"
boost_subset <- subset(boost, new_order %in% c("3", "4", "5"))
subset <- SubsetData(boost_subset, ident.use = c("Gaba", "Glut"))
genes <- c("CASP9", "CASP3", "CASP6")
p <- VlnPlot(subset, features = genes, pt.size = 0, cols = cols.cluster, combine = F)  #+NoLegend()+xlab('')#+theme(axis.text.x = element_text(angle = 0,  hjust = 0.5))
# patchwork <- wrap_plots(p, byrow = T)

3.14.1 CASP3

p[[2]] + # scale_y_continuous(breaks=c(0.0, 0.3,0.6, 0.9))+
theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20), 
    plot.title = element_text(size = 18)) + xlab("") + ylab("") + NoLegend()

3.14.2 CASP9

p[[1]] + scale_y_continuous(breaks = c(0, 0.3, 0.6, 0.9)) + theme(axis.text.x = element_text(angle = 0, 
    hjust = 0.5), axis.text = element_text(size = 20), plot.title = element_text(size = 18)) + 
    xlab("") + ylab("") + NoLegend()
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.

3.14.3 CASP6

p[[3]] + # scale_y_continuous(breaks=c(0.0, 1.0, 2.0, 3.0), labels=c('0.0', '1.0', '2.0',
# '3.0'))+
theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20), 
    plot.title = element_text(size = 18)) + xlab("") + ylab("") + NoLegend()

genes <- c("BIRC2", "XIAP", "AREL1", "SIVA1", "USP11", "BAX", "BCL2", "BAK1", "CYCS")
p <- VlnPlot(subset, features = genes, pt.size = 0, cols = cols.cluster, combine = F)  #+NoLegend()+xlab('')#+theme(axis.text.x = element_text(angle = 0,  hjust = 0.5))
# patchwork <- wrap_plots(p, byrow = T)

3.14.4 BIRC2

p[[1]] + scale_y_continuous(breaks = c(0, 0.3, 0.6, 0.9)) + theme(axis.text.x = element_text(angle = 0, 
    hjust = 0.5), axis.text = element_text(size = 20), plot.title = element_text(size = 18)) + 
    xlab("") + ylab("") + NoLegend()
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.

3.14.5 XIAP

p[[2]] + # scale_y_continuous(breaks=c(0.0, 0.3,0.6, 0.9))+
theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20), 
    plot.title = element_text(size = 18)) + xlab("") + ylab("") + NoLegend()

3.14.6 AREL1

p[[3]] + scale_y_continuous(breaks = c(0, 1, 2, 3), labels = c("0.0", "1.0", "2.0", 
    "3.0")) + theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20), 
    plot.title = element_text(size = 18)) + xlab("") + ylab("") + NoLegend()
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.

3.14.7 BCL2

p[[7]] + # scale_y_continuous(breaks=c(0.0, 0.3,0.6, 0.9))+
theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20), 
    plot.title = element_text(size = 18)) + xlab("") + ylab("") + NoLegend()

3.14.8 BAK1

p[[8]] + scale_y_continuous(breaks = c(0, 1, 2, 3), labels = c("0.0", "1.0", "2.0", 
    "3.0")) + theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20), 
    plot.title = element_text(size = 18)) + xlab("") + ylab("") + NoLegend()
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.

3.14.9 BAX

p[[6]] + scale_y_continuous(breaks = c(0, 0.3, 0.6, 0.9)) + theme(axis.text.x = element_text(angle = 0, 
    hjust = 0.5), axis.text = element_text(size = 20), plot.title = element_text(size = 18)) + 
    xlab("") + ylab("") + NoLegend()
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.

3.14.10 USP11

p[[5]] + # scale_y_continuous(breaks=c(0.0, 1.0, 2.0, 3.0), labels=c('0.0', '1.0', '2.0',
# '3.0'))+
theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20), 
    plot.title = element_text(size = 18)) + xlab("") + ylab("") + NoLegend()

3.14.11 SIVA1

p[[4]] + # scale_y_continuous(breaks=c(0.0, 1.0, 2.0, 3.0), labels=c('0.0', '1.0', '2.0',
# '3.0'))+
theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20), 
    plot.title = element_text(size = 18)) + xlab("") + ylab("") + NoLegend()

3.14.12 CYCS

p[[9]] + # scale_y_continuous(breaks=c(0.0, 1.0, 2.0, 3.0), labels=c('0.0', '1.0', '2.0',
# '3.0'))+
theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20), 
    plot.title = element_text(size = 18)) + xlab("") + ylab("") + NoLegend()

4 SessionInfo

sessionInfo()
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur 10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
##  [1] grid      parallel  stats4    stats     graphics  grDevices utils    
##  [8] datasets  methods   base     
## 
## other attached packages:
##  [1] monocle3_0.2.3.0            SingleCellExperiment_1.12.0
##  [3] SummarizedExperiment_1.20.0 GenomicRanges_1.42.0       
##  [5] GenomeInfoDb_1.26.1         MatrixGenerics_1.2.0       
##  [7] matrixStats_0.57.0          clusterProfiler_3.18.0     
##  [9] patchwork_1.1.0             forcats_0.5.0              
## [11] stringr_1.4.0               dplyr_1.0.2                
## [13] purrr_0.3.4                 readr_1.4.0                
## [15] tidyr_1.1.2                 tibble_3.0.4               
## [17] ggplot2_3.3.2               tidyverse_1.3.0            
## [19] org.Hs.eg.db_3.12.0         AnnotationDbi_1.52.0       
## [21] IRanges_2.24.0              S4Vectors_0.28.0           
## [23] Biobase_2.50.0              BiocGenerics_0.36.0        
## [25] viridis_0.5.1               viridisLite_0.3.0          
## [27] Seurat_3.2.2               
## 
## loaded via a namespace (and not attached):
##   [1] reticulate_1.18        tidyselect_1.1.0       RSQLite_2.2.1         
##   [4] htmlwidgets_1.5.2      BiocParallel_1.24.1    Rtsne_0.15            
##   [7] scatterpie_0.1.5       munsell_0.5.0          codetools_0.2-18      
##  [10] ica_1.0-2              future_1.20.1          miniUI_0.1.1.1        
##  [13] withr_2.3.0            colorspace_2.0-0       GOSemSim_2.16.1       
##  [16] knitr_1.30             rstudioapi_0.13        ROCR_1.0-11           
##  [19] tensor_1.5             DOSE_3.16.0            listenv_0.8.0         
##  [22] labeling_0.4.2         GenomeInfoDbData_1.2.4 polyclip_1.10-0       
##  [25] bit64_4.0.5            farver_2.0.3           downloader_0.4        
##  [28] parallelly_1.21.0      vctrs_0.3.5            generics_0.1.0        
##  [31] xfun_0.31              R6_2.5.0               graphlayouts_0.7.1    
##  [34] rsvd_1.0.3             DelayedArray_0.16.0    bitops_1.0-6          
##  [37] spatstat.utils_1.17-0  fgsea_1.16.0           assertthat_0.2.1      
##  [40] promises_1.1.1         scales_1.1.1           ggraph_2.0.4          
##  [43] enrichplot_1.10.1      gtable_0.3.0           globals_0.14.0        
##  [46] goftest_1.2-2          tidygraph_1.2.0        rlang_1.0.3           
##  [49] splines_4.0.3          lazyeval_0.2.2         broom_0.7.2           
##  [52] BiocManager_1.30.10    yaml_2.2.1             reshape2_1.4.4        
##  [55] abind_1.4-5            modelr_0.1.8           backports_1.2.0       
##  [58] httpuv_1.5.4           qvalue_2.22.0          tools_4.0.3           
##  [61] bookdown_0.21          ellipsis_0.3.1         jquerylib_0.1.4       
##  [64] RColorBrewer_1.1-2     ggridges_0.5.2         Rcpp_1.0.5            
##  [67] plyr_1.8.6             zlibbioc_1.36.0        RCurl_1.98-1.2        
##  [70] rpart_4.1-15           deldir_0.2-3           pbapply_1.4-3         
##  [73] cowplot_1.1.0          zoo_1.8-8              haven_2.3.1           
##  [76] ggrepel_0.8.2          cluster_2.1.0          fs_1.5.0              
##  [79] magrittr_2.0.1         data.table_1.13.2      DO.db_2.9             
##  [82] lmtest_0.9-38          reprex_0.3.0           RANN_2.6.1            
##  [85] fitdistrplus_1.1-1     hms_0.5.3              mime_0.9              
##  [88] evaluate_0.14          xtable_1.8-4           readxl_1.3.1          
##  [91] gridExtra_2.3          compiler_4.0.3         KernSmooth_2.23-18    
##  [94] crayon_1.3.4           shadowtext_0.0.7       htmltools_0.5.2       
##  [97] mgcv_1.8-33            later_1.1.0.1          speedglm_0.3-2        
## [100] lubridate_1.7.9.2      DBI_1.1.0              tweenr_1.0.1          
## [103] formatR_1.7            dbplyr_2.0.0           MASS_7.3-53           
## [106] Matrix_1.2-18          cli_2.2.0              igraph_1.2.6          
## [109] pkgconfig_2.0.3        rvcheck_0.1.8          plotly_4.9.2.1        
## [112] xml2_1.3.2             bslib_0.3.1            XVector_0.30.0        
## [115] rvest_0.3.6            digest_0.6.27          sctransform_0.3.1     
## [118] RcppAnnoy_0.0.17       spatstat.data_1.5-2    rmarkdown_2.14        
## [121] cellranger_1.1.0       leiden_0.3.5           fastmatch_1.1-0       
## [124] uwot_0.1.9             shiny_1.5.0            lifecycle_0.2.0       
## [127] nlme_3.1-150           jsonlite_1.7.1         fansi_0.4.1           
## [130] pillar_1.4.7           lattice_0.20-41        fastmap_1.1.0         
## [133] httr_1.4.2             survival_3.2-7         GO.db_3.12.1          
## [136] glue_1.4.2             spatstat_1.64-1        png_0.1-7             
## [139] bit_4.0.4              ggforce_0.3.2          stringi_1.5.3         
## [142] sass_0.4.1             blob_1.2.1             memoise_1.1.0         
## [145] irlba_2.3.3            future.apply_1.6.0